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>High  resolution  power  spectra  of  the  ZUrich  sunspot  numbers  for  the  period 
1849  to  1978  are  obtained  using  the  maximum  entropy  technique  of  Burg. 

The  monthly  means  are  first  smoothed  using  a  digital,  least-squares, 
band-pass  filter  with  193  weights,  and  then  decimated  to  obtain  a  more 
manageable  series.  An  overall  power  spectrum,  which  displays  multiple 
structure  and  harmonic  series  is  then  obtained.  In  order  to  explain  the 
complications  in  the  spectrum,  a  dynamic  spectrum  is  displayed.  This 
is  obtained  by  finding  the  spectra  or  data  samples  68  years  long, 
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^repeatedly  slipped  by  2  years  for  a  total  of  78  spectra.  The  spectra  are 
then  much  simpler  and  vary  smoothly  with  epoch.  The  nominal  11 -year 
line  varies  in  period  from  about  10  to  12f  years.  In  the  process  of* 
obtaining  a  power  spectrum  a  prediction  error  filter  is  derived.  This 
linear  filter  may  then  be  used  to  make  predictions  as  follows:  using  sets 
of  data  containing  5  solar  cycles  of  unsmoothed  monthly  values  beginning 
with  cycles  9  through  IS,  predictions  of  the  next  12  months  of  the  time 
series  are  made  and  compared  to  the  observed  values.  BMS  errors  vary 
between  5  and  33  and  lead  to  an  expectation  that  the  error  of  prediction 
for  cycle  21  will  be  about  20.  The  entire  cycle  21  is  then  predicted  using 
50,  100,  ISO  prediction  error  coefficients.  The  maximum  of  cycle  21  is 
predicted  to  be  130  ±  20  at  1980. 1  ±  0. 2. 


-  -V  ti\ 


Unclassified 


■BCUNITV  CL MStFIC ATtOtt  OP  TMI*  PA*tfW*«*l  OaM  SMWWO 


AFGL-TR-79-0104 


INTERNATIONAL  I ' r « • ; •  r i  1 1 r  Nu.  1 03 

SOLAR-TERRESTRIAL  PREDICTIONS 

PROCEEDINGS  AND  WORKSHOP  PROGRAM  l).u,-:  Pec  2?tJ97B 
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Paul  E.  Fougere 

Air  Force  Geophysics  Laboratory 
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High  resolution  power  spectra  of  the  Zurich  sunspot  numbers 
for  the  period  181)9  to  1978  are  obtained  using  the  maximum 
entropy  technique  of  Burg.  The  monthly  means  are  first  smoothed 
using  a  digital,  least-squares ,  band-pass  filter  with  193  weights, 
and  then  decimated  to  obtain  a  more  manageable  series.  An  over¬ 
all  power  spectrum,  which  displays  multiple  structure  and 
harmonic  series  is  then  obtained. 

In  order  to  explain  the  complications  in  the  spectrum,  a 
dynamic  spectrum  is  displayed.  This  is  obtained  by  finding  the 
spectra  of  data  samples  66  years  long,  repeatedly  slipped  by  2 
years  for  a  total  of  78  spectra.  The  spectra  are  then  much  simpler 
and  vary  smoothly  with  epoch.  The  nominal  11-year  line  varies 
in  period  from  about  10  to  12}  years. 

In  the  process  of  obtaining  a  power  spectrum  a  prediction 
error  filter  is  derived.  This  linear  filter  may  then  be  used  to 
make  predictions  as  follows;  using  sets  of  data  containing  9 
solar  cycles  of  unsmoothed  monthly  values  beginning  with  cycles 
9  through  15,  predictions  of  the  next  12  months  of  the  time  series 
are  made  and  compared  to  the  observed  values.  RMS  errors  vary 
between  5  and  33  and  lead  to  an  expectation  that  the  error  of 
prediction  for  cycle  21  will  be  about  20.  The  entire  cycle  21 
is  then  predicted  using  50,  100,  150  prediction  error  coeffi¬ 
cients.  The  maximum  of  cycle  21  is  predicted  to  be  130  ±  20 
at  1980.1  ±  0.2. 


Introduct ion 

The  subject  of  this  paper  is  the  time  series  of  Zurich  sunspot  numbers, 
R  .  The  basic  data  set  consists  of  monthly  mean  values  of  R  beginning 
in  1749  and  running  without  gaps  to  October  1978,  the  latest  available  value. 
See  Chernosky  and  Hagan  (1958)  for  the  earlier  data.  There  is  some 
smoothing  built  into  the  way  in  which  sunspot  counts  are  obtained,  because 
each  daily  count  uses  the  entire  visible  hemisphere,  which  is  rotating  with 
a  period  of  27  days,  or  a  rate  of  13  1/3  degrees  per  day.  Thus  the  daily 
count  is  really  a  131-day  running  mean  of  those  spots  and  groups  which 
could  be  observed  on  a  13  1/3  degree  tune  centered  at  central  medidian. 
Despite  the  low-pass  filter,  one  of  the  striking  features  of  the  series  of 
monthly  means  of  R^  is  the  large  amplitude,  high  frequency  noise  which  rides 
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FIGURE  1.  SOLAR  CYCLES  1-7.  DIGITAL  LOW-PASS  FILTER  WITH  13  WEIGHTS 
HAS  BEEN  USED.  EACH  CYCLE  BEGINS  AT  TIME  OF  SMOOTHED 
MINIMUM. 


on  top  of  the  other  striking  feature,  the  11-year  period.  Accordingly,  it 
has  become  customary,  for  example  in  the  monthly  Solar  Geophysical  Data 
bulletin,  to  further  smooth  the  monthly  averages  by  using  a  running  mean 
of  13  monthly  values  wi th  the  2  extreme  months  weighted  at  2. 

So  dominant  is  the  appearance  of  the  11-year  cycle  that  the  individual 
cycles  have  been  assigned  numbers  beginning  with  cycle  1  in  17^9  to  cycle 
21,  the  current  cycle,  wMch^star ted  in  mid-19/6.  Since  the  time  of  minimum 
can  be  located  with  r  accuracy,  the  cycles  begin  and  end  at  a 

minimum.  Figure  1  shows  the  first  7  cycles  with  the  time  origin  set  at  the 
time  of  minimum.  These  data  have  been  smoothed  using  a  13-point  digital 
least-squares,  low-pass  filter  following  the  method  of  Behannon  and  Ness 
(1966).  This  simple  filter  has  a  more  nearly  accurate  low-pass  response 
than  that  of  the  13“point  running  mean. 

Despite  the  use  of  the  filter,  the  extreme  variability  of  the  data  is 
evident.  All  of  the  cycles  arc  bumpy;  cycle  length,  from  minimum  to 
minimum,  varies  between  9  and  15  years;  time  from  minimum  to  maximum  varies 
from  about  2}  to  6  years,  and  the  size  of  the  maximum  varies  by  a  factor 
of  1*.  These  early  7  cycles  have  been  criticized  in  the  classical  paper  on 
sun  spot  prediction  by  McNish  and  Lincoln  ( 1 9f*9) .  According  to  that  paper 
the  early  data  are  considerably  less  reliable  and  even  belong  to  a  different 
statistical  population  from  the  modern  data  beginning  in  183^. 

Cycles  8  through  20,  the  modern  data,  are  shown  in  Figure  2.  For  these 
data,  the  length  of  the  cycle  varies  less,  between  10  and  12J  years.  The 


FIGURE  2.  SOLAR  CYCLES  8-20.  DIGITAL  LOW-PASS  FILTER  WITH  13  WEIGHTS 

HAS  BEEN  USED.  EACH  CYCLE  BEGINS  AT  TIME  OF  SMOOTHED  MINIMUM. 

time  to  maximum  is  now  about  2\  to  **J>  years,  but  t  hr-  height  of  the  maximum 
still  varies  by  a  factor  of  A:  from  about  50  to  more  than  200. 

The  extreme  variability  of  these  data  is  further  illustrated  in  Figure 
3,  showing  a  3~dimens ional  representation  of  all  20  cycles.  Cycle  19,  the 
largest  ever,  which  began  in  195**,  is  so  large  that  it  completely  hides 
cycle  20.  --  . 

Power  Spectra 

Now  we  would  like  to  determine  the  power  spectrum  of  this  time  series 
using  the  maximum  entropy  method  (Burg,  1975)  which  works  admirably  on 
rather  short-term  series.  The  original  2758  monthly  means  constitute  a 
rather  long  series.  We  can  easily  remedy  this  by  decimation  but  only  after 
smoothing  to  prevent  aliasing.  Accordingly,  a  low-pass  filter  with  a 
cutoff  period  of  2J  years  was  designed.  At  the  same  time,  a  high-pass 
filter  was  used  to  remove  the  mean  and  very  long  time  trends.  The  response 
of  the  resulting  band-pass  filter  is  shown  in  Figure  A.  The  digital  least- 
squares  filter  used  193  weights  to  achieve  this  response.  Notice  the 
filter  beginning  to  cut  off  at  about  2],  years,  as  designed,  so  that  the 
response  at  2  years  and  shorter  times  (high  frequencies)  is  close  to  zero. 
Thus  we  can  sample  the  output  safely  once  a  year  and  introduce  no  aliasing 
because  all  frequencies  higher  than  the  new  Nyguist  frequency  of  0.5  cycles 
per  year  have  been  removed. 


FIGURE  3. 


PERSPECTIVE  PLOT  OF  CYCLES  1-20. 


FIGURE  A.  RESPONSE  OF  DIGITAL  LEAST-SQUARES,  BAND-PASS  FILTER.  193 

WEIGHTS  ARE  USED.  FRCQUFNCY  IS  IN  UNITS  OF  CYCLES  PER  YEAR. 


The  results  are  shown  in  Figure  5.  The  original  unsmoothed,  monthly 
means  are  given  in  the  bottom  panel  with  values  running  from  0  to  250  and 
the  dates  running  from  the  year  1 7 SO  to  the  year  2000.  The  band-pass 
filtered  output  is  shown  in  the  top  panel,  and  the  differences  between  the 
original  and  band-pass  data,  the  residuals,  are  given  in  the  middle  panel. 
Thus  the  original  data  set  in  the  bottom  panel  is  the  sum  of  the  top  2 
panels,  month  by  month.  The  top  panel  shows  a  smooth,  zero-mean  time  series 
which  is  still  not  stationary  because  the  amplitudes  are  quite  variable. 

The  residual  series,  in  the  central  panel,  contains  essentially  all  of  the 
high-frequency  noise  plus  the  slowly-varying  mean  or  DC  level.  It  is  now 
quite  safe  to  decimate  the  band-pass  filtered  series  by  12  to  obtain  211 
yearly  values. 

The  maximum  entropy  method  of  Burg  (1975)  was  applied  to  these  211 
numbers  and  the  number  of  prediction  error  filter  weights,  analogous  to 
the  number  of  lags  in  a  Blackman-Tukey  spectrum,  was  varied  from  4  to  130. 
The  3"d imens ional  representation  of  these  spectra  is  shown  in  Figure  6. 

The  nominal  11-year  line,  which  carries  most  of  the  power,  splits  into  two 
very  stable  lines  with  periods  of  10.9  and  9-93  years.  There  is  also  quite 
clear  evidence  of  power  peaking  at  about  12  years  and  again  at  about  83 
years.  The  second  harmonic  of  the  11 -year  doublet  is  also  a  doublet  at 
about  5.5  and  4.8  years.  The  richness  of  structure  in  the  spectrum  is 
simply  another  manifestation  of  the  extreme  variability  of  the  time  series. 

The  complicated  spectrum  also  can  be  explained  by  taking  small  over¬ 
lapping  segments  of  data,  performing  spectral  analysis  on  each  and  watching 
the  spectrum  change  with  time.  Such  a  dynamical  spectrum  is  shown  in 
Figure  7.  Each  spectrum  is  based  on  a  segment  of  data  6>6  years  long. 
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FIGURE  5.  SUNSPOT  NUMBERS  VERSUS  DATA  FROM  17 50  TO  PRESENT  (OCT.  1978). 
TOP  PANEL:  BAND-PASS  FILTERED  DATA  (193  WEIGHTS). 

MIDDLE  PANEL:  ORIGINAL  MINUS  BAND-PASS. 

BOTTOM  PANEL:  ORIGINAL,  UNSM00THED  MONTHLY  MEANS. 


FREQUENCY  IN  CPY 

FIGURE  6.  PERSPECTIVE  PLOT  OF  MAXIMUM  ENTROPY  SPECTRA  OF  BAND-PASS  FILTERED 
DATA  SAMPLED  ONCE  PER  YEAR.  THE  NUMBER  OF  PREDICTION  ERROR 
FILTER  (PEF)  WEIGHTS  RANGES  BETWEEN  k  AND  130  IN  STEPS  OF  2. 
PERIODS  IN  YEARS  ARE  GIVEN  AT  APPROPRIATE  PEAKS. 


FREQUENCY  IN  CPY 

FIGURE  7.  DYNAMIC  MAXIMUM  ENTROPY  SPECTRUM  OF  66  YEARS  (APPROXIMATELY  6 
CYCLES)  OF  BAND-PASS  FILTERED  DATA  SAMPLED  ONCE  PER  YEAR. 

THE  START  YEAR  ADVANCES  BY  2  YEARS  FOR  EACH  SPECTRUM  SHOWN. 
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The  segments  are  overlapped  by  64  points,  so  that  the  start  year  moves  by 
2  years  for  each  new  segment.  We  are  now  concentrating  on  a  narrow  band 
of  periods  between  8  and  12J  years.  The  peaks  in  the  earliest  data  are  too 
small  to  be  seen  in  these  linear  power  spectral  density  plots.  The  nominal 
11-year  line  shows  an  instantaneous  period  varying  between  10  and  12J  years 
while  the  power  spectral  density  varies  by  over  two  orders  of  magnitude. 


Pred ict ions 


At  the  heart  of  the  maximum  entropy  method  is  the  determination  of  a 
prediction  error  filter.  This  is  a  linear  filter  whose  scaler  product 
with  a  segment  of  data  yields  the  error  in  a  one-step-ahead  prediction. 

We  find  the  filter  essentially  by  minimizing  the  sum  of  squares  of  such 
prediction  errors  over  the  set  of  observations.  Once  we  have  such  a  filter 
we  can  use  it  to  make  predictions  off  both  ends  of  the  series  as  far  as 
we  wish  to  go  simply  by  applying  the  filter  and  moving  over  the  data  and 
the  newly  acquired  predictions. 

The  predictions  of  the  sunspot  numbers  were  tested  as  follows.  Begin¬ 
ning  in  1843,  5  cycles  of  original  unsmoothed  monthly  values,  containing 
708  monthly  values,  constitute  a  data  Set  for  which  spectra  are  obtained 
for  a  number  of  prediction-error  filter-weights  ranging  from  SO  to  300  in 
steps  of  50.  For  each  filter,  predictions  were  made  for  the  following  12 
months  and  the  RMS  prediction  error  was  determined.  This  is  possible 
because  predictions  and  observed  values  which  hod  not  been  used  to  obtain 
the  predictions  are  available.  The  entire  procedure  was  repeated  seven 
times,  using  5  cycles  of  monthly  data  beginning  with  cycles  9,  10,  11 
through  cycle  15,  and  making  predictions  into  cycle's  14  through  20  respect¬ 
ively. 

The  mean  square  prediction  errors,  all  quite  reasonable  numbers  from 
a  low  of  5  to  a  high  of  33,  are  shown  in  Figure  8.  Here  are  plotted  the 
predicted  cycle  numbers  from  14  through  21  and  the  RMS  prediction  errors 
for  12  predictions  running  from  0  through  40.  Curiously,  the  RMS  error 
seems  to  alternate  in  magnitude,  being  low  for  even  cycles,  and  appreciably 
higher  for  odd-numbered  cycles.  Accordingly  a  guess  is  made  that  the  RMS 
prediction  error  for  the  real  predictions  for  cycle  21  will  be  around  20. 

Predictions  were  made  using  50  to  300  weights  for  the  entire  cycle 
21  and  plotted.  The  most  consistent  behavior  for  the  entire  cycle  comes 
from  sets  with  100,  150,  and  200  weights  and  these  are  shown  in  Figure  9- 
Smoothing  these  results,  the  predictions  of  the  maximum  and  its  date 
become : 


MAX  (R^)  ^  130  +  20 
DATE  =  1980.1  i  0.2 

The  entire  cycle  seems  quite  reasonable  with  a  minimum  around  1986.2. 

In  summary,  I  have  shown  some  pictures  of  the  sunspot  cycles  and 
discussed  their  extreme  variability.  I  have  filtered  these  data  with  a 
digital  band-pass  filter,  and  decimated  the  resulting  series.  1  have 
shown  maximum  entropy  power  spectra  varying  in  a  stable  manner  as  the 
number  of  weights  is  varied  and  also  a  dynamical  spectrum  showing  that  the 
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FIGURE  8.  RMS  ERRORS,  FOR  100,  150  AND  200  WEIGHTS.  FIVE  SOLAR  CYCLES  OF 

ORIGINAL  UNSMOOTHEO  MONTHLY  MEAN  ARE  USED  AS  INPUT  TO  THE  MAXIMUM 
ENTROPY  POWER  SPECTRAL  ANALYSIS  PROGRAM.  NUMBER  OF  OBSERVATIONS 
WAS  708,  708,  684,  684,  660.  648  AND  624  RESPECTIVELY. 


FIGURE  9.  PREDICTION  OF  UNSMOOTHEO  SUNSPOT  NUMBERS  MADE  USING  THREE 

DIFFERENT  SETS  OF  PREDICTION  ERROR  FILTERS  WITH  100,  150  AND 
200  WEIGHTS.  DATA  USED  ARE  670  ORIGINAL  UNSMOOTHED  MONTHLY 
MEANS  STARTING  IN  JANUARY,  1923. 


nominal  11-year  line  varies  in  period  from  10  to  12]  years.  I  have  tested 
a  prediction  technique  with  seven  separate  data  sets,  each  5  cycles  long, 
to  find  an  expected  RMS  prediction  error  of  about  20.  Finally  I  have 
made  predictions  for  the  remainder  of  cycle  21  which  should  have  a  maximum 
of  130  +  20  at  1980. 1  ♦  0.2. 
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